% Numerical example: K-12 broadband
% Authors: Gaurab Aryal, Pallavi Pal, Charlie Murry, and Arnab Palit
%
% Objective: Demonstrate how bundling reduces prices
% Focus: Two projects/schools with exogenous entry

%% Clear workspace and set random seed
clc
clear
rng(666)

%% Initialize parameters
s = 2;      % Number of schools
L = 1000;    % Number of simulations

%% Determine bids
Bids = []; 
Gamma = [4, 8, 12, 16]; 
ns     = [2, 3, 5, 10]; 

% Run simulations for different alphas and number of bidders
for gamma_id = 1:numel(Gamma)
    for n_id = 1:numel(ns)
        Bids = [BIDDING(ns(n_id), s, L, Gamma(gamma_id)); Bids];
    end
end

% save 
save('Bids.mat', 'Bids', '-v7.3'); 
%% Figures 
% Plot all bids and winning bids only for n=2 and create
% tables of averages across all 16 cases 
close all

% Define marker properties for each gamma value
marker_styles = {'o', 's', 'd', '^'};  % circle, square, diamond, triangle
marker_sizes = [50, 50, 50, 50];       % marker sizes
marker_colors = [0 0 0; 0.3 0.3 0.3; 0.5 0.5 0.5; 0.7 0.7 0.7]; % black to light gray
marker_faces = {[0 0 0], 'none', [0.5 0.5 0.5], [0.7 0.7 0.7]}; % filled/unfilled pattern

% Subsample parameter for Figure 1 (to reduce clutter)
subsample_fraction = 0.25;  % Use 25% of points for clearer visualization

% Figure 1 (a): All Total Bids 
fig1 = figure('Position', [100, 100, 900, 700]);
hold on

n = 2;
rng(42);  % Set seed for reproducible subsampling

for i = 1:4
    aa = 16 - (i-1) * 4;
    
    % Extract and prepare data
    bids_pre = Bids(aa).bids_pre;
    total_bids_pre = sum(bids_pre(:, 1:2, :), 2);
    total_bids_post = Bids(aa).bids_post;
    
    temp_pre = reshape(total_bids_pre, n*L, 1);
    temp_post = reshape(total_bids_post, n*L, 1);
    
    % Subsample for clearer visualization
    n_points = length(temp_pre);
    n_subsample = round(n_points * subsample_fraction);
    subsample_idx = randperm(n_points, n_subsample);
    
    temp_pre_sub = temp_pre(subsample_idx);
    temp_post_sub = temp_post(subsample_idx);
    
    % Plot with distinct markers
    scatter(temp_pre_sub, temp_post_sub, ...
        marker_sizes(i), ...
        marker_styles{i}, ...
        'MarkerEdgeColor', marker_colors(i,:), ...
        'MarkerFaceColor', marker_faces{i}, ...
        'LineWidth', 1.2, ...
        'MarkerEdgeAlpha', 0.7, ...
        'MarkerFaceAlpha', 0.6, ...
        'DisplayName', sprintf('\\Gamma = %d', Gamma(i)));
end

% Plot 45° reference line
plot([0 200], [0 200], 'k-', 'LineWidth', 2.5, ...
    'DisplayName', '45° line', 'HandleVisibility', 'off');

% Labels and formatting
xlabel('Bids Pre Bundling (b_1+b_2)', 'FontSize', 20, 'FontWeight', 'bold')
ylabel('Bids Post Bundling', 'FontSize', 20, 'FontWeight', 'bold')
title('All Total Bids [n=2]', 'FontSize', 22, 'FontWeight', 'bold')

% Grid and axes
grid on
set(gca, 'FontSize', 18, 'LineWidth', 1.5, 'GridAlpha', 0.3, ...
    'GridLineStyle', ':', 'Layer', 'top')
xlim([45, 200])
ylim([45, 200])
axis square

% Improved legend
leg = legend('Location', 'northwest', 'FontSize', 16, 'Box', 'on');
leg.Title.String = 'Cost Savings';
leg.Title.FontWeight = 'bold';

% Add text annotation for 45° line
text(180, 175, '45° line', 'FontSize', 14, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'right', 'BackgroundColor', 'white', ...
    'EdgeColor', 'k', 'Margin', 2);

hold off

% Save with high quality
set(fig1, 'PaperPositionMode', 'auto', 'Color', 'white');
print(fig1, 'bundling_comparison_all_bids_n2', '-depsc', '-r600');
saveas(fig1, 'bundling_comparison_all_bids_n2.fig');
exportgraphics(fig1, 'bundling_comparison_all_bids_n2.pdf', 'ContentType', 'vector');
exportgraphics(fig1, 'bundling_comparison_all_bids_n2.png', 'Resolution', 300);

% Figure 1(b): Winning Bids Only (Payments)
fig2 = figure('Position', [100, 100, 900, 700]);
hold on

% First plot 45° line so it's behind data
plot([0 200], [0 200], 'k-', 'LineWidth', 2.5, ...
    'DisplayName', '45° line');

for i = 1:4
    aa = 16 - (i-1) * 4;
    
    % Extract winning bids
    bids_pre = Bids(aa).bids_pre;
    bids_post = Bids(aa).bids_post;
    
    win_bid_pre = min(bids_pre(:,1,:), [], 1) + min(bids_pre(:,2,:), [], 1);
    win_bid_post = min(bids_post(:,1,:), [], 1);
    
    temp_pre = reshape(win_bid_pre, L, 1);
    temp_post = reshape(win_bid_post, L, 1);
    
    % Sort for better visualization
    combined = [temp_pre, temp_post];
    combined = sortrows(combined, 1);
    
    % Plot with distinct markers
    scatter(combined(:,1), combined(:,2), ...
        marker_sizes(i) + 10, ...  % Slightly larger for winning bids
        marker_styles{i}, ...
        'MarkerEdgeColor', marker_colors(i,:), ...
        'MarkerFaceColor', marker_faces{i}, ...
        'LineWidth', 1.3, ...
        'MarkerEdgeAlpha', 0.8, ...
        'MarkerFaceAlpha', 0.7, ...
        'DisplayName', sprintf('\\Gamma = %d', Gamma(i)));
end

% Labels and formatting
xlabel('Payment Pre Bundling', 'FontSize', 20, 'FontWeight', 'bold')
ylabel('Payment Post Bundling', 'FontSize', 20, 'FontWeight', 'bold')
title('Winning Bids (n=2)', 'FontSize', 22, 'FontWeight', 'bold')

% Grid and axes
grid on
set(gca, 'FontSize', 18, 'LineWidth', 1.5, 'GridAlpha', 0.3, ...
    'GridLineStyle', ':', 'Layer', 'top')
xlim([0, 200])
ylim([0, 200])
axis square

% Improved legend
leg = legend('Location', 'northwest', 'FontSize', 16, 'Box', 'on');
leg.Title.String = 'Cost Savings';
leg.Title.FontWeight = 'bold';

% Add text annotation for 45° line
text(190, 185, '45° line', 'FontSize', 14, 'FontWeight', 'bold', ...
    'HorizontalAlignment', 'right', 'BackgroundColor', 'white', ...
    'EdgeColor', 'k', 'Margin', 2);

hold off

% Save with high quality
set(fig2, 'PaperPositionMode', 'auto', 'Color', 'white');
print(fig2, 'bundling_comparison_winning_bids_n2', '-depsc', '-r600');
saveas(fig2, 'bundling_comparison_winning_bids_n2.fig');
exportgraphics(fig2, 'bundling_comparison_winning_bids_n2.pdf', 'ContentType', 'vector');
exportgraphics(fig2, 'bundling_comparison_winning_bids_n2.png', 'Resolution', 300);


%% Calculate averages
[avg_all_bids, avg_winning_bids] = calculate_averages(Bids);

%% BIDDING Function
function Output = BIDDING(n, s, L, Gamma)
    % BIDDING Simulates bidding in procurement auctions with cost complementarities
    % and when costs are Exponential. Won't work for others.
    %
    % Inputs:
    %   n - Number of bidders (positive integer)
    %   s - Number of schools/items (positive integer, typically 2)
    %   L - Number of simulations (positive integer)
    %   alpha - Complementarity parameter (scalar in [0,1])
    %
    % Outputs:
    %   Output - Struct containing:
    %     bids_pre: Array of bids for individual procurements (n x s x L)
    %     bids_post: Array of bids for bundled procurements (n x L)
    
    % Parameter initialization
    f_c        = makedist("Exponential", "mu", 40);  % Cost distribution (Exponential with mean 40)
    f_c        = truncate(f_c, 10, 100);             % Truncate c>10
    F          = @(c) cdf(f_c, c);                   % Cumulative distribution function of costs
    total_cost = @(c1, c2) c1 + c2 - Gamma;          % Total cost function with complementarities
    costs      = random(f_c, [n, s, L]);             % Generate random costs for all bidders, schools, and simulations
    max_cost   = max(costs(:));                      % Maximum cost across all draws
    
    % Individual procurements
    % Bidding function: b = c + markup - B
    % where markup = (1/(1-F(c)))^(n-1) * integral from c to c_max of (1-F(t))^(n-1) dt
    % B = Gamma x (1-F_c(cost of other school))^{(n-1)}

    B = @(c2) Gamma .* (1 - F(c2)).^(n-1);  % Adjustment term
    
    integral_func = @(c1) integral(@(t)(1-F(t)).^(n-1), c1, max_cost);  % Integral part of markup
    beta_pre_func = @(c1, c2) c1 + (1./(1-F(c1))).^(n-1) .* integral_func(c1) - B(c2);
    
    bids_pre = zeros(n, s, L);
    for l = 1:L
        cost = costs(:, :, l);
        % Calculate bids for each school using vectorized operations
        bids_pre(:, 1, l) = arrayfun(@(i) beta_pre_func(cost(i, 1), cost(i, 2)), 1:n);
        bids_pre(:, 2, l) = arrayfun(@(i) beta_pre_func(cost(i, 2), cost(i, 1)), 1:n);
    end
    
    % Bundled procurements (only for Exponential) 
    % Suppose every bidder has to submit total bid for both schools. Then the
    % objective function (bid_total - Total cost) x Pr(bid_total is minimum)
    % so the bidding function should have similar structure as before...
    % except the right cost cdf is the cdf of total cost
    
    F_tc = @(c) convolution_cdf(c, F, 10, 100, Gamma); 
    total_costs = total_cost(costs(:, 1, :), costs(:, 2, :));  % Calculate total costs
    max_total_cost = max(total_costs(:));  % Maximum total cost
    
    % Integral part of markup for bundled procurement
    integral_func_tc = @(c) integral(@(t)(1-F_tc(t)).^(n-1), c, max_total_cost);
    
    % Bidding function for bundled procurement
    beta_post = @(c) c + (1 ./ (1 - F_tc(c))).^(n-1) .* integral_func_tc(c);
    
    % Calculate bids for bundled procurement
    bids_post = arrayfun(beta_post, total_costs);
    
    % Prepare output
    Output.bids_pre = bids_pre;
    Output.bids_post = bids_post;
end

%% Function to calculate average bids and average winning bids
function [avg_all_bids, avg_winning_bids] = calculate_averages(Bids)
    ns = [2, 3, 5, 10];
    avg_all_bids = zeros(4, 4);
    avg_winning_bids = zeros(4, 4);
    
    for i = 1:4
        for j = 1:4
            aa = 16 - (i-1) * 4 - (j-1);
            n = ns(i);
            bids_pre = Bids(aa).bids_pre;
            bids_post = Bids(aa).bids_post;
            
            total_bids_pre = sum(bids_pre(:,1:2,:), 2);
            
            avg_all_bids(i,j) = mean([reshape(total_bids_pre, [], 1); 

            reshape(bids_post, [], 1)]);
            
            win_bid_pre = min(bids_pre(:,1,:), [], 1) + min(bids_pre(:,2,:), [], 1);
            
            win_bid_post = min(bids_post(:,1,:), [], 1);
            
            avg_winning_bids(i,j) = mean([reshape(win_bid_pre, [], 1); 

            reshape(win_bid_post, [], 1)]);
        end
    end
end